home *** CD-ROM | disk | FTP | other *** search
/ Tech Arsenal 1 / Tech Arsenal (Arsenal Computer).ISO / tek-19 / iritsm3s.zip / GEOMVALS.C < prev    next >
C/C++ Source or Header  |  1992-01-25  |  10KB  |  277 lines
  1. /*****************************************************************************
  2. *   "Irit" - the 3d polygonal solid modeller.                     *
  3. *                                         *
  4. * Written by:  Gershon Elber                Ver 0.2, Mar. 1990   *
  5. ******************************************************************************
  6. *  Module to evaluate geometric properties of the objects such as area,      *
  7. * volume, number of polygons etc.                         *
  8. *****************************************************************************/
  9.  
  10. #include <stdio.h>
  11. #include <ctype.h>
  12. #include <math.h>
  13. #include <string.h>
  14. #include "program.h"
  15. #include "allocate.h"
  16. #include "convex.h"
  17. #include "geomat3d.h"
  18. #include "geomvals.h"
  19. #include "objects.h"
  20. #include "windows.h"
  21. #include "graphgen.h"
  22.  
  23. static RealType PolygonXYArea(VertexStruct *VHead);
  24. static RealType Polygon3VrtxXYArea(PointType Pt1, PointType Pt2, PointType Pt3);
  25.  
  26. /*****************************************************************************
  27. *  Routine to evaluate the Area of the given geom. object, in object unit.   *
  28. * Algorithm (for each polygon):                    V3         *
  29. * 1. Set Polygon Area to be zero.                   /\         *
  30. *    Make a copy of the original polygon             /      \          *
  31. *    and transform it to a XY parallel plane.           /        \V2         *
  32. *    Find the minimum Y value of the polygon           V4/         |         *
  33. *    in the XY plane                     \         |         *
  34. * 2. Let V(0) be the first vertex, V(n) the last one       \         |         *
  35. *    for i goes from 0 to n-1 add to Area the area of         \_______|       *
  36. *    below edge V(i), V(i+1):                     V0      V1      *
  37. *    PolygonArea += (V(i+1)x - V(i)x) * (V(i+1)y' - V(i)y') / 2             *
  38. *    where V(i)y' is V(i)y - MinY, where MinY is the polygon minimum Y value *
  39. * 3. Note that the result of step 2 is the area of the polygon itself.       *
  40. *    However it might be negative, so take the absolute result of step 2 and *
  41. *    add it to the global ObjectArea.                         *
  42. * Note step 2 is performed by another aux. routine below: PolygonXYArea.     *
  43. *****************************************************************************/
  44. double
  45. #ifdef __MSDOS__
  46. cdecl
  47. #endif /* __MSDOS__ */
  48. PolyObjectArea(ObjectStruct *PObj)
  49. {
  50.     RealType ObjectArea = 0.0;
  51.     MatrixType RotMat;
  52.     PolygonStruct *Pl;
  53.     VertexStruct *V, *VHead;
  54.  
  55.     if (!IS_POLY_OBJ(PObj))
  56.     FatalError("Geometric property requested on non polygonal object");
  57.  
  58.     if (IS_POLYLINE_OBJ(PObj)) {
  59.     WndwInputWindowPutStr("Warning: geometric object is polyline - zero area");
  60.     return 0.0;
  61.     }
  62.  
  63.     Pl = PObj -> U.Pl.P;
  64.     while (Pl != NULL) {
  65.     V = VHead = CopyVList(Pl -> V);      /* Dont transform original object. */
  66.     /* Create the trans. matrix to transform the polygon to XY parallel  */
  67.     GenRotateMatrix(RotMat, Pl -> Plane);               /* plane. */
  68.     do {
  69.         MatMultVecby4by4(V -> Pt, V -> Pt, RotMat);
  70.  
  71.         V = V -> Pnext;
  72.     }
  73.     while (V != NULL && V != VHead);
  74.  
  75.     ObjectArea += PolygonXYArea(VHead);
  76.  
  77.     MyFree((char *) VHead, ALLOC_VERTEX);      /* Free the vertices list. */
  78.  
  79.     Pl = Pl -> Pnext;
  80.     }
  81.  
  82.     return ObjectArea;
  83. }
  84.  
  85. /*****************************************************************************
  86. *  Routine to evaluate the area of the given polygon projection on the XY    *
  87. * plane. Note the polygon does not have to be on a XY parallel plane, as     *
  88. * only its XY projection is considered (Z is ignored). Returns the area of   *
  89. * its XY parallel projection.                             *
  90. *  See GeomObjectArea above for algorithm:                     *
  91. *****************************************************************************/
  92. static RealType PolygonXYArea(VertexStruct *VHead)
  93. {
  94.     RealType PolygonArea = 0.0, MinY;
  95.     VertexStruct *V = VHead, *Vnext;
  96.  
  97.     MinY = V -> Pt[1];
  98.     V = V -> Pnext;
  99.     while (V != VHead && V != NULL /* Should not happen! */) {
  100.     if (MinY > V -> Pt[1]) MinY = V -> Pt[1];
  101.  
  102.     V = V -> Pnext;
  103.     }
  104.  
  105.     Vnext = V -> Pnext;
  106.     MinY *= 2.0;          /* Instead of subtracting twice each time. */
  107.     do {
  108.     /* Evaluate area below edge V-Vnext relative to Y level MinY. Note   */
  109.     /* it can come out negative, but thats o.k. as the sum of all these  */
  110.     /* quadraliterals should be exactly the area (up to correct sign).   */
  111.     PolygonArea += (Vnext -> Pt[0] - V -> Pt[0]) *
  112.                 (Vnext -> Pt[1] + V -> Pt[1] - MinY) / 2.0;
  113.     V = Vnext;
  114.     Vnext = V -> Pnext;
  115.     }
  116.     while (V != VHead && V != NULL /* Should not happen! */);
  117.  
  118.     return ABS(PolygonArea);
  119. }
  120.  
  121. /*****************************************************************************
  122. *   Routine to evaluate the Volume of the given geom object, in object unit. *
  123. *   This routine has a side effect that all non-convex polygons will be      *
  124. * splitted to convex ones.                             *
  125. * Algorithm (for each polygon, and let ObjMinY be the minimum OBJECT Y):     *
  126. *                                V3         *
  127. * 1. Set Polygon Area to be zero.                   /\         *
  128. *    Let V(0) be the first vertex, V(n) the last.         /      \          *
  129. *    For i goes from 1 to n-1 form triangles           /        \V2         *
  130. *    by V(0), V(i), V(i+1). For each such           V4/         |         *
  131. *    triangle do:                     \         |         *
  132. *    1.1. Find the vertex (out of V(0), V(i), V(i+1))      \         |         *
  133. *         with the minimum Z - TriMinY.                 \_______|       *
  134. *    1.2. The volume below V(0), V(i), V(i+1) triangle,         V0      V1      *
  135. *      relative to ObjMinZ level, is the sum of:                 *
  136. *      1.2.1. volume of V'(0), V'(i), V'(i+1) - the                 *
  137. *         area of projection of V(0), V(i), V(i+1) on XY parallel     *
  138. *         plane, times (TriMinZ - ObjMinZ).                 *
  139. *      1.2.2. Assume V(0) is the one with the PolyMinZ. Let V"(i) and     *
  140. *         V"(i+1) be the projections of V(i) and V(i+1) on the plane  *
  141. *         Z = PolyZMin. The volume above 1.2.1. and below the polygon *
  142. *         (triangle!) will be: the area of quadraliteral V(i), V(i+1),*
  143. *         V"(i+1), V"(i), times distance of V(0) for quadraliteral    *
  144. *         plane divided by 3.                         *
  145. *    1.3. If Z component of polygon normal is negative add 1.2. result to    *
  146. *      ObjectVolume, else subtract it.                     *
  147. *****************************************************************************/
  148. double
  149. #ifdef __MSDOS__
  150. cdecl
  151. #endif /* __MSDOS__ */
  152. PolyObjectVolume(ObjectStruct *PObj)
  153. {
  154.     int PlaneExists;
  155.     RealType ObjVolume = 0.0, ObjMinZ, TriMinZ, Area, PolygonVolume, Dist;
  156.     PointType Pt1;
  157.     PlaneType Plane;
  158.     PolygonStruct *Pl;
  159.     VertexStruct *V, *VHead, *Vnext, *Vtemp;
  160.  
  161.     if (!IS_POLY_OBJ(PObj))
  162.     FatalError("Geometric property requested on non polygonal object");
  163.  
  164.     if (IS_POLYLINE_OBJ(PObj)) {
  165.     WndwInputWindowPutStr("Warning: geometric object is polyline - zero area");
  166.     return 0.0;
  167.     }
  168.  
  169.     ObjMinZ = INFINITY;     /* Find Object minimum Z value (used as min level). */
  170.     Pl = PObj -> U.Pl.P;
  171.     while (Pl != NULL) {
  172.     V = VHead = Pl -> V;
  173.     do {
  174.         if (V -> Pt[2] < ObjMinZ) ObjMinZ = V -> Pt[2];
  175.         V = V -> Pnext;
  176.     }
  177.     while (V != VHead && V != NULL);
  178.     Pl = Pl -> Pnext;
  179.     }
  180.  
  181.     ConvexPolyObject(PObj);           /* Make sure all polygons are convex. */
  182.     Pl = PObj -> U.Pl.P;
  183.     while (Pl != NULL) {
  184.     PolygonVolume = 0.0; /* Volume below poly relative to ObjMinZ level. */
  185.     V = Vtemp = VHead = Pl -> V;/* We set VHead to be vertex with min Z: */
  186.     do {
  187.         if (V -> Pt[2] < Vtemp -> Pt[2]) Vtemp = V;
  188.         V = V -> Pnext;
  189.     }
  190.     while (V != VHead && V != NULL);
  191.     VHead = Vtemp;       /* Now VHead is the one with lowest Z in polygon! */
  192.     TriMinZ = VHead -> Pt[2];         /* Save this Z for fast access. */
  193.  
  194.     V = VHead -> Pnext;
  195.     Vnext = V -> Pnext;
  196.     do {
  197.         /* VHead, V, Vnext form the triangle - find volume 1.2.1. above: */
  198.         Area = Polygon3VrtxXYArea(VHead -> Pt, V -> Pt, Vnext -> Pt);
  199.         PolygonVolume += Area * (TriMinZ - ObjMinZ);
  200.  
  201.         /* VHead, V, Vnext form the triangle - find volume 1.2.2. above: */
  202.         Area = sqrt(SQR(V -> Pt[0] - Vnext -> Pt[0]) +   /* XY distance. */
  203.             SQR(V -> Pt[1] - Vnext -> Pt[1])) *
  204.            ((V -> Pt[2] + Vnext -> Pt[2]) / 2.0 - TriMinZ);
  205.         PT_COPY(Pt1, V -> Pt);
  206.         Pt1[2] = TriMinZ;
  207.         if ((PlaneExists =
  208.          CGPlaneFrom3Points(Plane, V -> Pt, Vnext -> Pt, Pt1)) == 0) {
  209.         /* Try second pt projected to Z = TriMinZ plane as third pt. */
  210.         PT_COPY(Pt1, Vnext -> Pt);
  211.         Pt1[2] = TriMinZ;
  212.         PlaneExists =
  213.             CGPlaneFrom3Points(Plane, V -> Pt, Vnext -> Pt, Pt1);
  214.         }
  215.         if (PlaneExists) {
  216.         Dist = CGDistPointPlane(VHead -> Pt, Plane);
  217.         PolygonVolume += Area * ABS(Dist) / 3.0;
  218.         }
  219.  
  220.         V = Vnext;
  221.         Vnext = V -> Pnext;
  222.     }
  223.     while (Vnext != VHead);
  224.  
  225.     if (Pl -> Plane[2] < 0.0)
  226.         ObjVolume += PolygonVolume;
  227.     else
  228.         ObjVolume -= PolygonVolume;
  229.  
  230.     Pl = Pl -> Pnext;
  231.     }
  232.  
  233.     return ObjVolume;
  234. }
  235.  
  236. /*****************************************************************************
  237. *  Routine to evaluate the area of the given triangle projected to the XY    *
  238. * plane, given as 3 Points. Only the X & Y components are considered.         *
  239. *  See PolyObjectArea above for algorithm:                     *
  240. *****************************************************************************/
  241. static RealType Polygon3VrtxXYArea(PointType Pt1, PointType Pt2, PointType Pt3)
  242. {
  243.     RealType PolygonArea = 0.0, MinY;
  244.  
  245.     MinY = MIN(Pt1[1], MIN(Pt2[1], Pt3[1])) * 2.0;
  246.  
  247.     PolygonArea += (Pt2[0] - Pt1[0]) * (Pt2[1] + Pt1[1] - MinY) / 2.0;
  248.     PolygonArea += (Pt3[0] - Pt2[0]) * (Pt3[1] + Pt2[1] - MinY) / 2.0;
  249.     PolygonArea += (Pt1[0] - Pt3[0]) * (Pt1[1] + Pt3[1] - MinY) / 2.0;
  250.  
  251.     return ABS(PolygonArea);
  252. }
  253.  
  254. /*****************************************************************************
  255. *  Routine to count number of polygons in given geometric object.         *
  256. *****************************************************************************/
  257. double
  258. #ifdef __MSDOS__
  259. cdecl
  260. #endif /* __MSDOS__ */
  261. PolyCountPolys(ObjectStruct *PObj)
  262. {
  263.     int Count = 0;
  264.     PolygonStruct *Pl;
  265.  
  266.     if (!IS_POLY_OBJ(PObj))
  267.     FatalError("Geometric property requested on non polygonal object");
  268.  
  269.     Pl = PObj -> U.Pl.P;
  270.     while (Pl != NULL) {
  271.     Count++;
  272.     Pl = Pl -> Pnext;
  273.     }
  274.  
  275.     return ((double) Count);
  276. }
  277.